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We have studied spin structures of fluctuation-driven fractionalized vortices and topological spin 
order in 2D nematic superfiuids of cold sodium atoms. Our Monte Carlo simulations suggest a 
softened 7r-spin disclination structure in a half-quantum vortex when spin correlations are short 
ranged; in addition, calculations indicate that a unique non-local topological spin order emerges 
simultaneously as cold atoms become a superfluid below a critical temperature. We have also 
estimated fluctuation-dependent critical frequencies for half-quantum vortex nucleation in rotating 
optical traps and discussed probing these excitations in experiments. 



Quantum number fractionalization has been one of the 
most fundamental and exciting concepts studied in mod- 
ern many-body physics and topological field theories [J 
0,11 11- During the past few years, low dimensional frac- 
tionalized quantum states have further been proposed 
to be promising candidates for carrying out fault tol- 
erant quantum computationjSjand their realizations in 
optical lattices were explored]! 0, 0] . A closely related 
topic in which there has also been a growing interest is 
vortex fractionalization in cold gases (see for instance 
Ref . d [H UH EH ) • Especially, motivated by experiments 
on low dimensional cold gases (l3j. Mukerjee et al stud- 
ied 2D superfiuids of cold atoms and analyzed the role 
played by fractionalized vortices in phase transitions fiol]. 
However, spin structures of those half-quantum vortices 
induced by thermal fluctuations and potential topolog- 
ical order 1A\ haven't been thoroughly explored and re- 
main to be understood. In this Letter we illustrate spin 
structures of fractionalized vortices; in addition, we also 
show that 2D quantum gases with short ranged spin cor- 
relations can have a topological spin order. We further 
estimate critical nucleation frequencies of fractionalized 
vortices in optical traps which can potentially be studied 
in experiments 15]. 

Our simulations illustrate that in a fundamental vortex 
carrying one- half of circulation quantum h/m (h is the 
Planck constant and m is the mass of atoms), there exists 
a softened spin disclination (i.e. a disclination in the ab- 
sence of spin stiffness) even when local spin moments are 
strongly fluctuating at finite temperatures. The topolog- 
ical winding number associated with softened spin discli- 
nations is conserved as far as the phase rigidity remains fi- 
nite. This effectively leads to a non-local topological spin 
order. Such a nonlocal order is absent in a conventional 
condensate of atoms or pairs of atoms. We have further 
studied creation of these excitations in rotating superfiu- 
ids and obtained fluctuation-dependent critical frequen- 
cies for half-quantum vortex(HQV) nucleation. HQVs in 
traps can be probed by measuring a precession of eige- 
naxes of surface quadrupole modes. 

We employ the Hamiltonian introduced previously for 



F=l sodium atoms in optical lattices (9l. [iH]. 
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Here k is the lattice site index and < kl > are the nearest 
neighbor sites. \x is the chemical potential and tr, is the 
one-particle hopping amplitude. Two coupling constants 
are b L (c L ) = b(c)^- J dr((f>* w (r)<j) w {r)) 2 ; b, c are effec- 
tive s-wave scattering lengths, <p w is the localized Wan- 
nier function for atoms in a periodical potential. Oper- 
ators f/'i, ce = x,y,z create hyperfine spin-one atoms in 
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respectively. The spin and number operators are defined 
as S a = —ie a ^t/)tipry , and p = vb\vb a . Spin correlations 

are mainly induced by interaction ClS 2 . Here we consider 
antiferromagnetic spin-dependent interactions such as in 
sodium atoms where cl > 0. Minimization of this antifer- 
romagnetic spin-dependent interaction requires that the 
order parameter ^ a (—< vb\. >) be a real vector up to a 

global phase, i.e. ^ = y/Nnexp(ix) where n is a unit di- 
rector on a two-sphere, exp(i\) represents a phase direc- 
tor and N is the number of atoms per site. All low energy 
degrees of freedom are characterized by configurations 
where n and \ vary slowly in space and time[9j. Low 
lying collective modes include spin-wave excitations with 
energy dispersion u)(q) = v s q, v s — ^fcjj^a and phase- 
wave excitations with w(<?) = v p q, v p — \fbjjyia (here a 
is the lattice constant). In one-dimensions, low energy 
quantum fluctuations destroy spin order leading to quan- 
tum spin disordered superfiuids [13] • In two-dimensions, 
the amplitude of quantum spin fluctuations is of order of 
cl/^l and is negligible in shallow lattices as is order-of- 
magnitude bigger than cl- At finite temperatures, spin 
correlations are mainly driven by long wave length ther- 
mal fluctuations, analogous to quantum ID cases. This 
aspect was also paid attention to previously and normal- 
superfluid transitions were investigated [To| . 

We therefore study the following Hamiltonian which 
effectively captures long wave length thermal fluctuations 
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FIG. 1: (color online) a)Spin(circle) and phase(square) cor- 
relation lengths versus temperatures; inset is for the renor- 
malized phase coupling constant J p (in units of J). Dashed 
lines are the fits to exp(A/ \JT — T c ) for phase correla- 
tions with T c m 0.35 J and exp(B/T) for spin correlations, 
b) Thermal average (W a ) (circle) and (W p ) (square) versus 
the loop perimeter at T = 0.33J(solid line) and T = 
0.45J(dashed line). c)-f): (W s ) hv (circle), (W p }hv (square) 
and (C) hv (diamond) averaged over configurations with HQV 
boundary conditions. c)-d),f) are for T/J = 0.45,0.33,0.20 re- 
spectively; in e), we also show {W a , P )hv normalized in terms of 
back ground values (W s ,p}bg at T = 0.33 J. (W s ,p)bg (dashed 
fines) are averaged over configurations with uniform phases 
at the boundary. 
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here states at each site are specified by two unit directors: 
a nematic director, n = (cos <po sin 9q, sin 4>o sin 0q, cos 0q) 
and a phase director, 3? = (cos x, sin x). Jki = 2NtL 
is the effective coupling between two neighboring sites 
and depends on N, the number of atoms per site. The 
model is invariant under the following local Ising gauge 
transformation: rij — ► s^rii, <fr — > Sj$j, and Si = ±1. 
In the following, we present results of our simulations on 
2D superfluids, especially spin structures, energetics of 
HQVs and nucleation of HQVs in rotating optical traps 
using the effective Hamiltonian in Eq[21 



HQV and underlying topological spin order Around a 
HQ V, both phase and nematic directors rotate slowly by 
180°; in polar coordinates (0,p), a HQV in condensates 
is represented by ^(9,p) = ^jN(p) exp(i#/2)n(0), with 
n = (cos(0/2),sin(0/2),O)[H[lZl. The question here is 
whether, when nematic directors are not ordered, a spin 
disclination is still present in a HQV. To fully take into 
account 2D thermal fluctuations, we carry out Monte 
Carlo simulations on a square lattice of 128 x 128 sites 
and study spatial correlations between a HQ V and a tt- 
spin disclination, and topological order. 

We first identify critical temperatures of the normal- 
superfluid phase transition by calculating correlations 
and the phase rigidity. The gauge-invariant quadrupolc- 
quadrupole correlation functions we have studied are 
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/ s 'nri,r 2 )=<Q^(r 1 )Q^(r 2 )>. 

Here Q s a p{?\) = ni, a m,^ - (l/3)6 a p, a = x,y,z; 

Qaj3( r i) = *i,a*i,/3 - (1/ 2 )<W a = x,y. In simula- 
tions, we have studied these correlation functions and 
found that the phase correlation length for / p (ri, be- 
comes divergent at a temperature 0.35 J which is identi- 
fied as a critical temperature T c . We also calculate the 
phase rigidity or the renormalized phase coupling J p 
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here S\ is a small phase difference applied across the op- 
posite boundaries of the lattice and F is the correspond- 
ing free energy. We indeed find that it approaches zero 
at T c while at T = J p takes a bare value J. Mean- 
while, the spin correlation function / s (ri,r2) remains to 
be short ranged across T c . By extrapolating our data 
to lower temperatures, we find that the spin correlation 
length diverges only at T — (see FigfT^,). Our simula- 
tions for correlation lengths are in agreement with previ- 
ous results in Ref . [Toj ; they are also consistent with the 
continuum limit of the model in Eqf2] which is equivalent 
to an XY model and an 0(3) nonlinear-sigma model. 

In order to keep track of the winding of nematic direc- 
tors in a wildly fluctuating back ground, we introduce the 
following gauge invariant 7r-rotation checking operator, 
which is essentially a product of sign-checking operators 
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Here the product is carried out along a closed square- 
shape path C centered at the origin of a 2D lattice. W SiP 
can be either +1 or —1; and W^Wp) is —1 when C en- 
closes a 7r-spin disclination (HQV). The gauge invariant 
circulation of supercurrent velocity (in units of irh/m) 
is defined as C = \ J^ec sign(n k ■ nj) sin(xfc - Xi)\ this 
quantity is equal to one in a HQ V. 
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In our simulations, we investigate the winding number 
< Ws,p >hv averaged over configurations where phase 
directors rotate by 180° around the boundary of the lat- 
tice and the center plaquette. At temperatures above 
the normal-supcrfluid transition temperature T c , both 
winding numbers W s . p and circulation C are averaged 
to zero within our numerical accuracy (see FigfT]). And 
our choice of boundary conditions does not lead to a vor- 
tex or disclination configuration in the absence of phase 
rigidity. Below T c) the circulation C is averaged to one in- 
dicating that the boundary conditions effectively project 
out HQV configurations. Meanwhile, we observe loop- 
perimeter dependent (W s>p )hv which can be attributed 
to the background fluctuations of HQ V or disclination 
pairs. The loop-perimeter dependence of (W s , P )hv here is 
almost identical to that for uniform boundary conditions, 
i.e. the back ground value. After normalizing {W s p }hv 
in terms of background winding numbers {W SjP )t,g, we 
find both and ffi^ approach -1 (see FigHJ. 

We thus demonstrate that a softened disclination is spa- 
tially correlated with a HQV. At the temperatures we 
carry out these simulations the spin correlation length is 
sufficiently short compared to the size of the lattice. At 
further lower temperatures, the spin correlation length 
becomes longer than the lattice size and fluctuations of 
pairs of disclination-anti-disclination are strongly sup- 
pressed; < W s ,p >hv are equal to —1 for almost all loops, 
which corresponds to a mean field result. 

Results in Fig[T] indicate that there exists a softened 
spin disclination in a HQV. This is a distinct feature 
in our systems and there exist no such additional mag- 
netic structures in HQ Vs in conventional molecular con- 
densates of atom pairs discussed perviouslyfUj]. Thus, 
7r-disclinations like HQ Vs have logarithmically divergent 
energies and are fully suppressed in ground states. Our 
results also illustrate that although the average local 
spin quadrupole moments Q s a p vanish because of strong 
fluctuations, an overall 7r-rotation of nematic directors 
in disclinations is still conserved because of a coupling 
to the superfluid component. This coupling between a 
HQV and disclination can also be attributed [19[ to a cou- 
pling between Higgs matter and discrete gauge fields [20l|. 
Furthermore, the absence of unbound 7r-disclinations in 
superfluids indicates a topological order, similar to the 
one introduced previously for an isotropic phase of liq- 
uid crystal 21] . Consequently, once a conventional phase 
order appears below a critical temperature, a topological 
spin order simultaneously emerges while spin correlations 
remain short ranged. 

The emergent topological order can be further verified 
by examining the average of product-operator W StP over 
all configurations (with open boundaries). Above the 
normal-superfluid transition temperature T c , we again 
find that W s>p both are averaged to zero within our 
numerical accuracy implying proliferation of unbound 
HQVs or disclinations. Below T c , we study the loop- 
perimeter dependence of average winding numbers (W s , p ) 
and find that both ln(W p ) and \n(W s ) are linear functions 



of loop-perimeter analogous to the Wilson-loop-product 
of deconfining gauge fields [IH; if there were unbound 
disclinations, one should expect that ln(W p ) is propor- 
tional to, instead of the loop-perimeter, the loop-area 
which represents the number of unbound disclinations 
enclosed by the loop. 

Critical frequency for HQV nucleation Let us now 
turn to the nucleation of those excitations in rotating 
traps [H [H [H, [H, M, [H, H3 . To understand the criti- 
cal frequency for nucleation, we study the free energy of a 
vortex, in a rotating frame, as a function of the distance 
r from the axis of a cylindrical optical trap (the axis is 
along the z-direction) , 

F h . v .(r)=Fl v Xr)-flL z (r). (6) 

Here F° v (r) is the free energy of a HQ V located at dis- 
tance r from the trap axis in the absence of rotation, 
L z (r) is the angular momentum of the vortex state and 
£1 is the rotating frequency. 

In a 2D lattice without a trapping potential, F® v is 
approximately equal to \ (J p + J s )ln(L/a), with lead- 
ing contributions from phase winding and spin twisting; 
here J P;S are renormalized phase and spin coupling re- 
spectively and L is the size of system. For an integer- 
quantum vortex (IQV), F® is equal to irJ p \n(L/a). The 
ratio between and F® depends on the ratio J s / J p or 
spin fluctuations; in the limit L approaches infinity, the 
ratio F® v /F% changes discontinuously from i at T = 
where J p « J s — J to \ at finite low temperatures in 2D 
where J s vanishes. In simulations of a finite trap (see be- 
low) , because of a finite size effect we find that this ratio 
varies from 0.5 to 0.2 smoothly as temperatures increase 
from to T c . 

To study nucleation of half-quantum vortices in an op- 
tical trap, we assume a nearly harmonic trapping po- 
tential V(r) = l/2mu 2 r r 2 , with ui tr being the trap fre- 
quency. The average number of particles per site N(r) 
has a Thomas-Fermi profile; N(r) — N (l — r 2 /R^ F ), 
here Nq is the number density at the center and Rtf is 
the Thomas-Fermi radius. Furthermore, the optical lat- 
tice potential along the axial direction is sufficiently deep 
so that atoms are confined in a two-dimensional xy plane; 
the in-plane lattice potential depth is set to be 5Er(Er 
is the photon recoil energy). 

For the trap and lattice geometry described above, we 
calculate parameters in Eq[T]and obtain ti, = 77nK, hi, — 
187nK and cl = lOnK. For Nq — 1.2 and trap frequency 
Lo tr — 120Hz, we also find that Rtf = 43a and Lh.o. = 
4o where Lh.o.{ = 1 / ^J2mu)tr) 1S the harmonic oscillator 
length. The coupling Jki in Eq[2]depends on the distance 
from the center of trap and at the center, the coupling 
is about 154nK. In non-rotating or slowly rotating traps, 
the free energy maximum is located at the center and 
there should be no vortices in the trap. As frequencies 
are increased, a local energy minimum appears at the 
center and becomes degenerate with the no-vortex state 
at a thermodynamic critical frequency (which is about 
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FIG. 2: (color online) a) The free energy of a HQV at a 
distance r from the trap center at rotation frequencies 0, 0.08, 
0.23 and 0.32 (in units of trap frequency cutr = 120 Hz) at T = 
0.33J. Details near the edge are shown in the inset. At the 
center, the exchange coupling is J = 154 nK. b) The critical 
frequencies f2 c for HQ Vs (solid line) and IQ Vs (dashed line) . 
Inset is the temperature-dependence of the ratio between the 
HQV energy F°.„, and the IQV energy F°. 

0.08uJtr at T = 0.33J); however because of a large energy 
barrier separating the two degenerate states as shown in 
Fig|2j vortices are still prohibited from entering the trap. 

Further speeding up rotations results in an energeti- 
cally lower and spatially narrower barrier. Within the 
range of temperatures studied, thermal activation turns 
out to be insignificant within an experimental time scale 
(~ 100ms) because of low attempt frequencies. So only 
when the spatial width of barrier becomes comparable to 
a hydrodynamic breakdown length [30l|. the barrier can 
no longer be felted and vortices start to penetrate into 
the trap. We use this criterion to numerically determine 
the dynamical critical frequency for vortex nucleation 
f2,P; for IQVs, the calculated f2^ is a flat function of T 
(see Figl2b) which is qualitatively consistent with earlier 
estimates [33j. 

For HQ Vs, F$ depends on the renormalized spin cou- 
pling J s and therefore the amplitude of spin fluctuations. 
Because of this, fi,P varies from about 0.32o7t r at T = 
where J s « J p and Fh. v . ~ 0.5F V due to a finite size 
effect (see the inset of Figj2]), to about 0.17<x>t r at tem- 
peratures close to T c where J s — and Fh. v . ~ 0.2F V . In 
other words, this unique temperature dependence can be 



considered to be an indicator of fluctuation-driven vor- 
tex fractionalization. It is worth remarking that in the 
thermodynamic limit where J s approaches zero at any 
finite temperatures, Fh. v , approaches jF v as mentioned 
before. Consequently, (about 0.17cj tr ) for HQVs is 
about one-half of the critical frequency for IQVs (about 
0.36wt r for the finite trap studied here). Also note that 
the zero temperature estimate of fif is close to the pre- 
viously obtained value of critical frequencies of HQ Vs in 
Bose-Einstein condensates 

The interaction between two HQ Vs with the same vor- 
ticity at a separation distance d contains two parts. One, 
V cc {> 0) is from interactions between two supercurrent 
velocity fields which is logarithmic as a function of d; 
and the other, V ss is from interaction between two spin 
twisting fields accompanying HQVs. For a disclination- 
anti disclination pair, in the dilute limit one finds that 
V cc ~ — V ss resulting in a cancellation of long range in- 
teractions. The resultant short-range repulsions lead to 
square vortex lattices found in numerical simulations [32J- 
For fluctuation-driven fractionalized vortices, V ss is al- 
most zero and the overall interactions are always loga- 
rithmically repulsive. HQVs nucleated in a rotating trap 
should therefore form a usual triangular vortex lattice. 

Individual vortex lines can be probed either by study- 
ing a precession of eigenaxes of surface quadrupole mode 
in rotating superfluids [33l] . In the later approach, one 
studies the angular momentum carried per particle in a 
HQ V state. When a HQV is nucleated in the trap, su- 
perfluids are no longer irrotational and the angular mo- 
mentum per particle is h/2 rather than h per particle 
for an integer vortex state. When a surface quadrupole 
oscillation across a rotating superfluid is excited, larger 
axes of quadrupole oscillation start to precess just as in 
the case of integer vortices. However, the precession rate 
is only one half of the value for an integer vortex state 
which can be studied in experiments. 

In conclusion, 2D superfluids of sodium atoms have 
a non-local topological spin order. In rotating traps, 
fluctuation-driven fractionalized vortices can nucleate at 
a critical frequency which is about half of that for integer 
vortices. Observation of these exotic excitations could 
substantially improve our understanding of topological 
order and fractionalization. We thank J. Zhang and Z. 
C. Gu for contributions at an early stage of the project. 
This work is in part supported by the office of the Dean 
of Science, UBC, NSERC (Canada), Canadian Institute 
for Advanced Research, and the A. P. Sloan foundation. 
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